1 Set working dir

here::i_am("4_cell_type_annotation.Rmd")
here::set_here()
print(paste("Current working directory:", here::here()))
## [1] "Current working directory: C:/Users/flyku/Documents/GitHub/NOTCH1-scRNAseq/rmarkdown"

2 Load data

source("functions.R")
#dir.create('../sample_obj')
combined <- qs::qread('combined.qsave')

DefaultAssay(combined) <- "RNA"
Idents(combined) <- combined$orig.ident

custom_color <-
  as.character(palette36.colors(36)[-2])
two_color <- c('#C0C0C0', '#B00D23')

provided_marker <- read_csv("../marker/marker_by_day.csv") %>%
  mutate(cell_type = factor(cell_type))

provided_marker
DimPlot(
  combined,
  reduction = "umap",
  cols = custom_color,
  label = T,
  pt.size = 0.4,
  repel = T,
  label.box = T
)

3 Annotate by each sample

3.1 Con0

Con0: all iPSCs

# Settings
this_sample_name <- "Con0"
this_day <- "day0"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if(file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <- FindClusters(this_combined, resolution = 0.2)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 1,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

tmp_ident <- as.factor(this_combined$seurat_clusters)
levels(tmp_ident) <- rep("iPSCs", length(levels(tmp_ident)))
this_combined <- AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

VlnPlot(this_combined,
        features = unique(this_markers),
        ncol = 1)

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

table(this_combined$cell_type)
## 
## iPSCs 
##  2825
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
## iPSCs 
##     1

3.2 N1KO0

N1KO0: all iPSCs

# Settings
this_sample_name <- "N1KO0"
this_day <- "day0"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if (file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <- FindClusters(this_combined, resolution = 0.2)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 1,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

tmp_ident <- as.factor(this_combined$seurat_clusters)
levels(tmp_ident) <- rep("iPSCs", length(levels(tmp_ident)))
this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

VlnPlot(this_combined,
        features = unique(this_markers),
        ncol = 1)

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

table(this_combined$cell_type)
## 
## iPSCs 
##  1556
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
## iPSCs 
##     1

3.3 Con2

Con2: iPSCs, Mesoderm, Cardiac mesoderm

# Settings
this_sample_name <- "Con2"
this_day <- "day2"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if (file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
  
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <-
    FindClusters(this_combined, resolution = 2.6, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

#VlnPlot(this_combined, features = this_markers, ncol = 3)
tmp_ident <- as.factor(this_combined$seurat_clusters)
tmp_levels <- rep("Mesoderm", 19)
tmp_levels[c(8, 11, 14, 16, 19)] <- "iPSCs"
tmp_levels[c(7, 15)] <- "Cardiac mesoderm"

levels(tmp_ident) <- tmp_levels


this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

VlnPlot(this_combined,
        features = unique(this_markers),
        ncol = 3)

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

table(this_combined$cell_type)
## 
##         Mesoderm Cardiac mesoderm            iPSCs 
##             2154              274              624
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
##         Mesoderm Cardiac mesoderm            iPSCs 
##        0.7057667        0.0897772        0.2044561

3.4 N1KO2

N1KO2: iPSCs, Mesoderm, Cardiac mesoderm

# Settings
this_sample_name <- "N1KO2"
this_day <- "day2"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if (file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <-
    FindClusters(this_combined, resolution = 0.3, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

#VlnPlot(this_combined, features = this_markers, ncol = 3)
tmp_ident <- as.factor(this_combined$seurat_clusters)

tmp_levels <- rep("Mesoderm", 4)
tmp_levels[c(4)] <- "iPSCs"
tmp_levels[c(1, 3)] <- "Cardiac mesoderm"

levels(tmp_ident) <- tmp_levels

this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

VlnPlot(this_combined,
        features = unique(this_markers),
        ncol = 3)

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

table(this_combined$cell_type)
## 
## Cardiac mesoderm         Mesoderm            iPSCs 
##             1790              813              131
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
## Cardiac mesoderm         Mesoderm            iPSCs 
##       0.65471836       0.29736650       0.04791514

3.5 Con5

Con5: iPSCs, Mesoderm, Cardiac mesoderm, SHF progenitor, FHF progenitor, Epicardial progenitor

# Settings
this_sample_name <- "Con5"
this_day <- "day5"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if (file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <-
    FindClusters(this_combined, resolution = 1.2, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

deg_filename <- paste0("../sample_obj/", this_sample_name, "_marker.csv")
if (!file.exists(deg_filename)) {
  cts_markers <-
    FindAllMarkers(this_combined,
                   only.pos = T,
                   min.pct = 0.5,
                   logfc.threshold = 0.3) %>%
    filter(pct.1 > 0.25)
  write.csv(cts_markers,
            paste0("../sample_obj/", this_sample_name, "_marker.csv"))
}

tmp_ident <- as.factor(this_combined$seurat_clusters)
#length(levels(tmp_ident))
tmp_ident <- as.factor(this_combined$seurat_clusters)
tmp_levels <- rep("Undecided", length(levels(tmp_ident)))
levels(tmp_ident) <- tmp_levels

this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
## Undecided 
##         1
VlnPlot(this_combined, features = unique(this_markers), ncol = 3)

3.6 N1KO5

N1KO5: iPSCs, Mesoderm, Cardiac mesoderm, SHF progenitor, FHF progenitor, Epicardial progenitor

# Settings
this_sample_name <- "N1KO5"
this_day <- "day5"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if(file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <- FindClusters(this_combined, resolution = 2, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

#VlnPlot(this_combined, features = unique(this_markers), ncol = 3)
deg_filename <- paste0("../sample_obj/", this_sample_name, "_marker.csv")
if (!file.exists(deg_filename)) {
  cts_markers <-
    FindAllMarkers(this_combined,
                   only.pos = T,
                   min.pct = 0.5,
                   logfc.threshold = 0.3)
  write.csv(cts_markers,
            paste0("../sample_obj/", this_sample_name, "_marker.csv"))
}


tmp_ident <- as.factor(this_combined$seurat_clusters)
tmp_levels <- rep("Cardiac mesoderm", length(levels(tmp_ident)))
tmp_levels[c(14)] <- "iPSCs"
tmp_levels[c(13)] <- "SHF progenitor"

levels(tmp_ident) <- tmp_levels


this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
## Cardiac mesoderm   SHF progenitor            iPSCs 
##       0.94006515       0.03387622       0.02605863
VlnPlot(this_combined, features = unique(this_markers), ncol = 3)

3.7 Con10

Con10: Cardiac mesoderm, SHF progenitor, FHF progenitor, Epicardial progenitor, Atrial cardiomyocytes, Ventricular cardiomyocytes, Pacemaker cells, Vascular smooth muscle cells, Cardiac fibroblast

and early cardiomyocytes

# Settings
this_sample_name <- "Con10"
this_day <- "day10"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if(file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <- FindClusters(this_combined, resolution = 0.8, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

#VlnPlot(this_combined, features = unique(this_markers), ncol = 3)
deg_filename <- paste0("../sample_obj/", this_sample_name, "_marker.csv")
if (!file.exists(deg_filename)) {
  cts_markers <-
    FindAllMarkers(this_combined,
                   only.pos = T,
                   min.pct = 0.5,
                   logfc.threshold = 0.3)
  write.csv(cts_markers,
            paste0("../sample_obj/", this_sample_name, "_marker.csv"))
}

tmp_ident <- as.factor(this_combined$seurat_clusters)
tmp_levels <- rep("Early cardiomyocytes", length(levels(tmp_ident)))
tmp_levels[c(9, 11, 14)] <- "FHF progenitor"
tmp_levels[c(3, 5)] <- "SHF progenitor"
tmp_levels[c(6)] <- "Epicardial progenitor"

levels(tmp_ident) <- tmp_levels

this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

#Con10
table(this_combined$cell_type)
## 
##  Early cardiomyocytes        SHF progenitor Epicardial progenitor 
##                  1212                   394                   179 
##        FHF progenitor 
##                   263
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
##  Early cardiomyocytes        SHF progenitor Epicardial progenitor 
##            0.59179688            0.19238281            0.08740234 
##        FHF progenitor 
##            0.12841797
VlnPlot(this_combined, features = unique(this_markers), ncol = 3)

3.8 N1KO10

N1KO10: Cardiac mesoderm, SHF progenitor, FHF progenitor, Epicardial progenitor, Atrial cardiomyocytes, Ventricular cardiomyocytes, Pacemaker cells, Vascular smooth muscle cells, Cardiac fibroblast

and early cardiomyocytes

# Settings
this_sample_name <- "N1KO10"
this_day <- "day10"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if(file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <- FindClusters(this_combined, resolution = 1, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

#VlnPlot(this_combined, features = unique(this_markers), ncol = 3)
deg_filename <- paste0("../sample_obj/", this_sample_name, "_marker.csv")
if (!file.exists(deg_filename)) {
  cts_markers <-
    FindAllMarkers(this_combined,
                   only.pos = T,
                   min.pct = 0.5,
                   logfc.threshold = 0.3)
  write.csv(cts_markers,
            paste0("../sample_obj/", this_sample_name, "_marker.csv"))
}

tmp_ident <- as.factor(this_combined$seurat_clusters)
tmp_levels <- rep("Epicardial progenitor", length(levels(tmp_ident)))
tmp_levels[c(12, 13)] <- "FHF progenitor"
tmp_levels[c(2, 7, 8, 9)] <- "SHF progenitor"
tmp_levels[c(3, 5, 14)] <- "Early cardiomyocytes"

levels(tmp_ident) <- tmp_levels


this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

#N1KO10
table(this_combined$cell_type)
## 
## Epicardial progenitor        SHF progenitor  Early cardiomyocytes 
##                  1116                   883                   573 
##        FHF progenitor 
##                   139
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
## Epicardial progenitor        SHF progenitor  Early cardiomyocytes 
##            0.41165622            0.32571007            0.21136112 
##        FHF progenitor 
##            0.05127259
VlnPlot(this_combined, features = unique(this_markers), ncol = 3)

3.9 Con14

Con14: SHF progenitor, FHF progenitor, Epicardial progenitor, Atrial cardiomyocytes, Ventricular cardiomyocytes, Pacemaker cells, Vascular smooth muscle cells, Cardiac fibroblast

# Settings
this_sample_name <- "Con14"
this_day <- "day14"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if(file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <- FindClusters(this_combined, resolution = 1, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

#VlnPlot(this_combined, features = unique(this_markers), ncol = 3)
deg_filename <- paste0("../sample_obj/", this_sample_name, "_marker.csv")
if (!file.exists(deg_filename)) {
  cts_markers <-
    FindAllMarkers(this_combined,
                   only.pos = T,
                   min.pct = 0.5,
                   logfc.threshold = 0.3)
  write.csv(cts_markers,
            paste0("../sample_obj/", this_sample_name, "_marker.csv"))
}

tmp_ident <- as.factor(this_combined$seurat_clusters)
tmp_levels <- rep("Early cardiomyocytes", length(levels(tmp_ident)))
tmp_levels[c(15)] <- "Vascular smooth muscle cells"
tmp_levels[c(1, 7, 9, 12, 13)] <- "Cardiac fibroblast"
levels(tmp_ident) <- tmp_levels

this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")

Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

table(this_combined$cell_type)
## 
##           Cardiac fibroblast         Early cardiomyocytes 
##                          481                          920 
## Vascular smooth muscle cells 
##                           29
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
##           Cardiac fibroblast         Early cardiomyocytes 
##                   0.33636364                   0.64335664 
## Vascular smooth muscle cells 
##                   0.02027972
VlnPlot(this_combined, features = unique(this_markers), ncol = 3)

3.10 N1KO14

N1KO14: SHF progenitor, FHF progenitor, Epicardial progenitor, Atrial cardiomyocytes, Ventricular cardiomyocytes, Pacemaker cells, Vascular smooth muscle cells, Cardiac fibroblast

# Settings
this_sample_name <- "N1KO14"
this_day <- "day14"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if(file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <- FindClusters(this_combined, resolution = 1, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

#VlnPlot(this_combined, features = unique(this_markers), ncol = 3)
deg_filename <- paste0("../sample_obj/", this_sample_name, "_marker.csv")
if (!file.exists(deg_filename)) {
  cts_markers <-
    FindAllMarkers(this_combined,
                   only.pos = T,
                   min.pct = 0.5,
                   logfc.threshold = 0.3)
  write.csv(cts_markers,
            paste0("../sample_obj/", this_sample_name, "_marker.csv"))
}

tmp_ident <- as.factor(this_combined$seurat_clusters)
tmp_levels <- rep("Early cardiomyocytes", length(levels(tmp_ident)))
tmp_levels[c(10)] <- "Epicardial progenitors"
tmp_levels[c(5, 8, 12, 16)] <- "Cardiac fibroblast"
tmp_levels[c(11, 15)] <- "Vascular smooth muscle cells"
levels(tmp_ident) <- tmp_levels



this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

table(this_combined$cell_type)
## 
##         Early cardiomyocytes           Cardiac fibroblast 
##                         1802                          496 
##       Epicardial progenitors Vascular smooth muscle cells 
##                          139                          138
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
##         Early cardiomyocytes           Cardiac fibroblast 
##                   0.69980583                   0.19262136 
##       Epicardial progenitors Vascular smooth muscle cells 
##                   0.05398058                   0.05359223
VlnPlot(this_combined, features = unique(this_markers), ncol = 3)

3.11 Con30

Con30: Atrial cardiomyocytes, Ventricular cardiomyocytes, Pacemaker cells, Vascular smooth muscle cells, Cardiac fibroblast

# Settings
this_sample_name <- "Con30"
this_day <- "day30"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if(file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <- FindClusters(this_combined, resolution = 2.6, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

#VlnPlot(this_combined, features = unique(this_markers), ncol = 3)
deg_filename <- paste0("../sample_obj/", this_sample_name, "_marker.csv")
if (!file.exists(deg_filename)) {
  cts_markers <-
    FindAllMarkers(this_combined,
                   only.pos = T,
                   min.pct = 0.5,
                   logfc.threshold = 0.3)
  write.csv(cts_markers,
            paste0("../sample_obj/", this_sample_name, "_marker.csv"))
}

tmp_ident <- as.factor(this_combined$seurat_clusters)
tmp_levels <- rep("Ventricular cardiomyocytes", length(levels(tmp_ident)))
tmp_levels[c(4)] <- "Atrial cardiomyocytes"
tmp_levels[c(11)] <- "Pacemaker cell"
tmp_levels[c(12)] <- "Cardiac fibroblast"
tmp_levels[c(14, 16)] <- "Vascular smooth muscle cells"

levels(tmp_ident) <- tmp_levels


this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

table(this_combined$cell_type)
## 
##   Ventricular cardiomyocytes        Atrial cardiomyocytes 
##                         1183                          132 
##               Pacemaker cell           Cardiac fibroblast 
##                           64                           61 
## Vascular smooth muscle cells 
##                           99
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
##   Ventricular cardiomyocytes        Atrial cardiomyocytes 
##                   0.76868096                   0.08576998 
##               Pacemaker cell           Cardiac fibroblast 
##                   0.04158545                   0.03963613 
## Vascular smooth muscle cells 
##                   0.06432749
VlnPlot(this_combined, features = unique(this_markers), ncol = 3)

3.12 N1KO30

N1KO30: Atrial cardiomyocytes, Ventricular cardiomyocytes, Pacemaker cells, Vascular smooth muscle cells, Cardiac fibroblast

# Settings
this_sample_name <- "N1KO30"
this_day <- "day30"

this_markers <- provided_marker %>%
  filter(str_detect(day, this_day)) %>%
  pull(gene)

Idents(combined) <- combined$orig.ident

# Load pre-computed object if exist
obj_filename <- paste0("../sample_obj/", this_sample_name, ".qsave")
if (file.exists(obj_filename)) {
  this_combined <- qs::qread(obj_filename)
} else {
  this_combined <- subset(combined, idents = this_sample_name)
  this_combined <- FindVariableFeatures(this_combined, verbose = F)
  this_combined <- RunPCA(this_combined, verbose = F)
  this_combined <- RunUMAP(this_combined, dims = 1:20, verbose = F)
  this_combined <- FindNeighbors(this_combined, verbose = F)
  this_combined <-
    FindClusters(this_combined, resolution = 2, verbose = F)
}


Idents(this_combined) <- this_combined$seurat_clusters
DimPlot(this_combined, reduction = "umap", label = T, label.box = T)

FeaturePlot(
  this_combined,
  features = unique(this_markers),
  ncol = 3,
  keep.scale = NULL,
  pt.size = 0.6,
  cols = two_color
)

#VlnPlot(this_combined, features = unique(this_markers), ncol = 3)
deg_filename <- paste0("../sample_obj/", this_sample_name, "_marker.csv")
if (!file.exists(deg_filename)) {
  cts_markers <-
    FindAllMarkers(this_combined,
                   only.pos = T,
                   min.pct = 0.5,
                   logfc.threshold = 0.3)
  write.csv(cts_markers,
            paste0("../sample_obj/", this_sample_name, "_marker.csv"))
}

tmp_ident <- as.factor(this_combined$seurat_clusters)
tmp_levels <- rep("Atrial cardiomyocytes", length(levels(tmp_ident)))
tmp_levels[c(8, 12, 13, 14)] <- "Ventricular cardiomyocytes"
tmp_levels[c(11)] <- "Pacemaker cell"
tmp_levels[c(2, 6)] <- "Cardiac fibroblast"
tmp_levels[c(10)] <- "Vascular smooth muscle cells"

levels(tmp_ident) <- tmp_levels



this_combined <-
  AddMetaData(this_combined, tmp_ident, col.name = "cell_type")


Idents(this_combined) <- this_combined$cell_type
DimPlot(this_combined, reduction = "umap")

# Save
qs::qsave(this_combined, paste0("../sample_obj/", this_sample_name, ".qsave"))

table(this_combined$cell_type)
## 
##        Atrial cardiomyocytes           Cardiac fibroblast 
##                          818                          261 
##   Ventricular cardiomyocytes Vascular smooth muscle cells 
##                          241                           69 
##               Pacemaker cell 
##                           64
table(this_combined$cell_type) / sum(as.numeric(table(this_combined$cell_type)))
## 
##        Atrial cardiomyocytes           Cardiac fibroblast 
##                   0.56297316                   0.17962836 
##   Ventricular cardiomyocytes Vascular smooth muscle cells 
##                   0.16586373                   0.04748796 
##               Pacemaker cell 
##                   0.04404680
VlnPlot(this_combined, features = unique(this_markers), ncol = 3)

4 UMAP of cell types using all samples

obj1 <- qread("../sample_obj/Con0.qsave")
obj2 <- qread("../sample_obj/Con2.qsave")
obj3 <- qread("../sample_obj/Con5.qsave")
obj4 <- qread("../sample_obj/Con10.qsave")
obj5 <- qread("../sample_obj/Con14.qsave")
obj6 <- qread("../sample_obj/Con30.qsave")
obj7 <- qread("../sample_obj/N1KO0.qsave")
obj8 <- qread("../sample_obj/N1KO2.qsave")
obj9 <- qread("../sample_obj/N1KO5.qsave")
obj10 <- qread("../sample_obj/N1KO10.qsave")
obj11 <- qread("../sample_obj/N1KO14.qsave")
obj12 <- qread("../sample_obj/N1KO30.qsave")


obj <-
  merge(
    x = obj1,
    y = c(obj2,
          obj3,
          obj4,
          obj5,
          obj6,
          obj7,
          obj8,
          obj9,
          obj10,
          obj11,
          obj12)
  )

combined <- AddMetaData(combined, obj$cell_type, col.name = "cell_type")
Idents(combined) <- combined$cell_type

DimPlot(combined, reduction = "umap", cols = custom_color)
Idents(combined) <- combined$cell_type
DimPlot(combined,
        reduction = "umap",
        cols = custom_color,
        pt.size = 0.8)

qs::qsave(combined,'combined.qsave')